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Non-diagonalizable and non-divergent susceptibility tensor in the Hamiltonian 
mean-field model with asymmetric momentum distributions 
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We investigate response to an external magnetic field in the Hamiltonian mean-field model, which 
is a paradigmatic toy model of a ferromagnetic body and consists of plane rotators like the XY spins. 

Due to long-range interactions, the external field drives the system to a long-lasting quasistationary 
state before reaching thermal equilibrium, and the susceptibility tensor obtained in the quasista¬ 
tionary state is predicted by a linear response theory based on the Vlasov equation. For spatially 
homogeneous stable states, whose momentum distributions are asymmetric with zero-means, the 
theory reveals that the susceptibility tensor for an asymptotically constant external field is neither 
symmetric nor diagonalizable, and the predicted states are not stationary accordingly. Moreover, 
the tensor has no divergence even at the stability threshold. These theoretical findings are confirmed 
by direct numerical simulations of the Vlasov equation for the skew-normal distribution functions. 

PACS numbers: 05.20.Dd, 05.70.Jk, 74.25.N- 


I. INTRODUCTION 

Long-range Hamiltonian systems have many remark¬ 
able features PJ, and one of them is existence of quasis¬ 
tationary states (QSSs) in the way of relaxation to ther¬ 
mal equilibrium. The lifetime of QSSs diverges with the 
number of particles consisting of the system Si , and 
hence QSSs are solely observable in a system with large 
population like self-gravitating systems 0. Dynamics 
of such a system is described by the Vlasov equation, or 
the collisionless Boltzmann equation, in the limit of large 
population @- 0 , and the QSSs, including thermal equi¬ 
librium states, are regarded as stable stationary solutions 
to the Vlasov equation. The system slowly goes towards 
thermal equilibrium with large but finite population due 
to finite size effects 0. 

The QSSs are observed not only in isolated systems, 
but also in systems under external fields. The initial 
QSS, which may and may not be in thermal equilib¬ 
rium, is driven to another QSS by the external field, and 
the resulting QSS is not necessarily in thermal equilib¬ 
rium. As a result, response to the external field may dif¬ 
fer from one obtained by statistical mechanics. Indeed, 
in the ferromagnetic model so-called Hamiltonian mean- 
held (HMF) model [9, 1(|, the critical exponents are ob¬ 
tained as 7 _ = 1/4 11] and S = 3/2 fl2l ] with the aid 
of a linear EH EH and a nonlinear EH response theories 
based on the Vlasov description respectively, while sta¬ 
tistical mechanics gives y_ = 1 and 5 = 3. Interestingly, 
with another exponent /3 = 1/2, the non-classical expo¬ 
nents satisfy the classical scaling relation q_ = /3(<5 — 1), 
and have universality on initial reference families of QSSs 
in a wide class of ID mean-held models fl5j . 

The universality is derived under the assumption that 
the initial distribution functions depend on position and 
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momentum only through the one-particle Hamiltonian 
with referring to the Jeans theorem jlhl] . Thus, the ini¬ 
tial states are symmetric with respect to momentum. 
The symmetric initial states are also used in studies on 
nonequilibrium statistical mechanics [l7l - |20f . the core¬ 
halo description of QSSs [2l[, nonequilibrium dynam¬ 
ics [22], and correlation and diffusion [23]. See also 
Refs 42 [24 1. 

Nevertheless, asymmetric momentum distributions ap¬ 
pear in beam-plasma systems (see |2jl l2Sl for instance), 
and are experimentally created in an ultracold plasma by 
optical pumping [29| . In the HMF model, such distribu¬ 
tions are stationary even asymmetric, and it is, therefore, 
natural to ask the response in the asymmetric case for 
completing the response theory. The main purpose of this 
article is to investigate the linear response against asymp¬ 
totically constant external field around spatially homoge¬ 
neous but asymmetric distributions in the HMF model. 
It is worth noting that, despite its simpleness, the model 
shares similar dynamics with the free-electron laser [3Cj 
and an anisotropic Heisenberg model under classical spin 
dynamics jjV • 

The HMF model consists of plane rotators like the XY 
spins, and the susceptibility tensor in the HMF model is 
of size 2x2 corresponding to the x- and y-directions of the 
rotators. For symmetric homogeneous states, the suscep¬ 
tibility tensor is directly diagonalized and experiences a 
divergence at the critical point of the second order phase 
transition, which is dynamically interpreted as the sta¬ 
bility threshold of the homogeneous states EMI. We 
then ask the two questions for asymmetric momentum 
distributions with zero-means: Is the susceptibility ten¬ 
sor symmetric and diagonalizable ? Does the response 
diverge at the stability threshold ? We will answer these 
questions negatively. The non-diagonalizable response 
tensor implies that the external field for x-direction in¬ 
duces the magnetization for y-direction, and such a re¬ 
sponse is unavoidable even changing the coordinate. Due 
to this non-diagonalizability, the predicted stat is not sta- 
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tionary, while the constant external field may drive the 
system to a stationary state asymptotically. In other 
words, the non-diagonalizability provides an example of 
discrepancy between the asymptotic states by the lin¬ 
ear dynamics and the full Vlasov dynamics. The non¬ 
divergence of response suggests that 7 + = 0 and (5=1, 
and interestingly, the scaling relation 7 + = j3{5— 1) holds, 
although j3 might be not well defined since the spatially 
inhomogeneous stationary states must be symmetric by 
the Jeans theorem 53 ■ ' 

This article is organized as follows. The HMF model 
and the linear responses are reviewed in Sec|TTl As an 
example of a family of asymmetric distributions, we in¬ 
troduce the skew-normal distributions, and investigate 
their stability in Sec lIIIl Theoretical consequences are 
examined by direct numerical simulations of the Vlasov 
equation in SeclIVI We discuss on stationarity of the pre¬ 
dicted state in Sec|V] The last section [VI] is devoted to a 
summary and discussions. 


II. THE HAMILTONIAN MEAN-FIELD MODEL 
AND LINEAR RESPONSE THEORY 

A. The model 

The HMF model with the time-dependent external 
magnetic field h = (h x (t),h y (t)) is expressed by the 
Hamiltonian 

N p 2 1 N 

H N (q,p,t) = ^ y + yrr Hi 1- cos fe “ ?fc)] 


2 2 TV 

3 =1 j»«=1 

N 


(i) 


— ^^[h x (t) cos qj + h y (t ) sin^j]. 
i=i 


The corresponding one-particle Hamiltonian is defined on 
the p space, which is (—7r, 7r] x 1R, as 


B[f]{q,p,t) = y - (M x + h x ) cosq — (M y + h y )smq 

( 2 ) 

where the magnetization vector ( M x ,M y ) is defined by 

(M x ,M y ) = j J (cos sin q)f(q,p, t)dqdp. (3) 

The one-particle distribution function / is governed by 
the Vlasov equation 


d£ 

dt 


Wf]J} = 0 , 


with the Poisson bracket defined by 

Of d g df dg 




dp dq dq dp 


( 4 ) 


( 5 ) 


One can straightforwardly check that any spatially ho¬ 
mogeneous states, /o(p), are stationary if the external 
field h is absent. 

We prepare a homogeneous stable stationary state 
fo(p) for t < 0, and add a small external field h for t > 0. 
To avoid an artificial rotation, we require the zero-mean 
for fo(p), and consider an asymptotically constant exter¬ 
nal field accordingly. For instance, we set 

(‘;D= e< ‘>(t) (6) 

using the Heaviside step function ©(f), and the external 
field drives the initial state /o to / = /o + fi asymptot¬ 
ically. Accordingly, the one-particle Hamiltonian %[/] 
changes from H 0 to H 0 + H x , where 

p 2 

H 0 = Hi = —(M liX + h x )cosq- {M 1>y + h y ) sin q 

( 7 ) 

and 

M hx = (cos q) 1 , M XiV = (sin q) 1 . ( 8 ) 

We introduced the averages of an observable B with re¬ 
spect to /o and f x as 

(B)j = jj B(q,p)fj(q,p)dqdp, ( j = 0,1). (9) 


B. Isothermal linear response 


It might be instructive to review the isothermal linear 
response to compare it with the Vlasov linear response 
theory which will be presented in the next subsection 

m 

The thermal equilibrium states of the HMF model are 
describe by the one-particle distribution functions of 

g-^fZo+ih) 

= ff^e-^o+nddqdp (10) 

Hereafter /3 represents not one of the critical exponents 
mentioned in Sec|l] but the inverse temperature. Ex¬ 
panding / into the power series of H\ and picking up to 
the linear order, we have 


(B) 1 = -p[(BHi) 0 -(B) 0 (H 1 ) 0 }. (11) 

Substituting cos q and sin q into B, we have the matrix 
formula 



where the correlation matrix C = (C va ) (y,cr £ {x, y}) 
is defined by 


((cos q cos q) Q (cos q sin q) 0 \ 
Wsinqcosg ) 0 (singsing ) 0 J 


(13) 
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Thus, the formal solution is 

(«:;;) = 11 - (£;)■ <i4) 

and the susceptibility tensor % = (Xvcr) {v,<r £ {x,j/}) 
defined by M = in the limit ||h|| 0 is 

X = [1 - cr'c. (15) 

Divergence of x appears at the critical point satisfying 
det(l — C) = 0. 

It is easy to show that the correlation matrix is now 
expressed by C = (/3/2)/ 2 , where / 2 is the 2x2 unit 
matrix. The susceptibility tensor is hence diagonalized 
and the diagonal elements are 

.. _ 0/2 T c 

Xxx — Xyy ~ 1 _ pj2 ~ J 1 — J 1 

with the critical temperature T c = 1/2 of the second 
order phase transition 0- The vanishing off-diagonal 
elements come from spatial homogeneity of fo(p), and 
symmetry of fo(p) is not necessary. 


C. Vlasov linear response 


The nonlinear response theory [l2} includes the linear 
response theory [13|,ll4j for symmetric fg(p) and provides 
a simple expression of the linear response [151 ] , but asym¬ 
metric fg(p) is out of range. Thus, we revisit the linear 
response theory. 

We introduce the Laplace transform defined by 

/•OO 

u(ui) = / u(t)e lut dt. (17) 

Jo 


The linear response theory gives the Laplace transform 
of Mi tV (t)), denoted by (M ltX (u), M 1%y {u)), as 


( A/i,x(w)\ 
\Ml,y(u})J 




(18) 


where the elements of matrix F = (F ucr ) are 

Fxx M = ~ir [ (—— + —r~) fo(p) d P 
2 JL \P — w P + uJ 

F xy(u) = ~p~ [ (— - fo(p) d P ( 19 ) 

2 * VP- w P + u) 

Fyx (w) — F X y(lU ) 

Fyy (w) = F xx (u). 

See the Appendix]^] for derivations. The integral contour 
L is the real p axis for Im(w) > 0, but is continuously 
modified for Irri(w) < 0 to avoid t he p oles at p = ±w by 
following the Landau’s procedure 0] ■ 

Temporal evolution of Mi, y ) is determined by 

performing the inverse Laplace transform, which picks up 


singularities of its Laplace transform (fl8l) . For instance, a 
pole at uj l gives a term having exp(—which implies 
the Landau damping for Irn(wL) < 0. Assuming that 
the reference fo(p) is stable, we have no singularities on 
the upper half w plane. Existence of singularities on the 
real axis of w is accidental for [1 — A(w)] _1 .F(u;), and 
we omit it. Then, the main singularity comes from the 
Heaviside step function of the external field (Fiji, whose 
Laplace transform is 


/h x (w)\ = -1 fh x \ 

\h y (uj)) iw \hy) 


( 20 ) 


Asymptotic values of Mi tX and M^ y are, therefore, ob¬ 
tained by picking up the pole at ui = 0 [14[, and 


(M hx {t)\ 

\Ml,y(t)) 



(t Oo), 


( 21 ) 


where the susceptibility tensor x = (Xw) is written in a 
similar form with (fT5l) as 

X = [1 - F(0)]~ 1 F(0). (22) 


Let us rewrite the above Vlasov susceptibility x by 
using the dispersion function 

D{u) = 1 + 7T [ ^Md p, w £ C. (23) 
JlP-u 

In the following we consider real u> which gives 

D(u) = 1 + 7T PV f ^Md p + i7r 2 /oM, w e I, 

J —00 P ^ 

(24) 

where PV represents the principal value. The dispersion 
function rewrites the susceptibility as 

1 /Re(D(0))-|D(0 )| 2 -Im(J9(0)) ^ 

X “ |D(0 )| 2 ^ lm(£>(0)) Re(£>(0)) - \D(0)\ 2 J ' 

(25) 

When fo(p) is symmetric and hence /g(0) = 0, imply¬ 
ing Im(D(0)) = 0 accordingly, the susceptibility tensor x 
is diagonal, and the diagonal elements are 

1 -Z>( 0 ) 

Xxx = Xyy = (26) 

as reported in Refs. 00- The susceptibility, there¬ 
fore, diverges at the point D( 0) = 0 corresponding to 
the stability threshold 0,0. O 11 the other hand, when 
/o(0) 7 ^ 0, the imaginary part of D( 0) does not van¬ 
ish and hence the susceptibility tensor K5l) enjoys two 
interesting features: (i) The tensor is neither symmetric 
nor diagonalizable by the real coordinate transformation, 
since the eigenvalues are not real, (ii) No divergence ap¬ 
pears even at the stability threshold, since |D( 0)| 2 > 0 . 
We note that, for homogeneous symmetric distributions, 
D( 0 ) > 0 is the stability criterion and hence the di¬ 
vergence appears at the stability threshold. However, 
D( 0) >0 is no more the stability criterion for the asym¬ 
metric case. A stability criterion for the asymmetric case 
will be introduced in Sec lIII Bl 
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III. SKEW-NORMAL DISTRIBUTION AND 
STABILITY 

A. Skew-normal distribution 

We introduce the skew-normal distribution for exam¬ 
ining the linear response theory and confirming the two 
features mentioned in Sec lII Cl Advantages of the skew- 
normal distribution are that it has the single peak which 
makes the stability criterion simpler, and that the ana¬ 
lytically obtained mean value helps to set the total mo¬ 
mentum zero. 

The density of skew-normal distribution is defined by 



- 3 - 2-10123 

x 


/snOu A, ix, a) = ^</> ■ 


where 


and 


4>(ir) 



(27) 


(28) 


(29) 


The parameter A represents the skewness, and A = 0 
results to the normal distribution. The mean value is 


FIG. 1. (color online) Skew-normal distributions with zero 
means and a = 1. A = -2,-1, 0,1 and 2, whose maximum 
points are from right to left. /o(0) is positive (resp. negative) 
for negative (resp. positive) A. 


the Nyquist method (H, HI}, which was applied to asym¬ 
metric double-peak distributions in the HMF model [35| . 

In our setting, the Nyquist method provides the sta¬ 
bility criterion as 


fo(p; A) has an exponentially 


xfsNdx 


p + crS 



A 

bTTW' 


(30) 


We test the homogeneous stationary states of the form 
job; A, lx, a) = t^/sn(p; A, p,a), (31) 

which is normalized as JJ .fodqdp = 1. To set the total 
momentum zero, we put 



Hereafter we fix the parameter <r as <r = 1. Then, the 
unique free parameter is the skewness A, and the distri¬ 
bution is simply denoted by job; A). Let p = q be the 
unique extreme point (the maximum point) depending 
on A. Some examples of the skew-normal distribution 
functions are exhibited in Fig [I] 


B. Nyquist method of stability 

For symmetric distributions fn ( p ), the formal stability 
criterion has been established [3f as 

fo{p) is formally stable D(0) > 0, (33) 

where D is the dispersion function (l24l) . To obtain the 
formal stability, job) is assumed as a function of one- 
particle Hamiltonian, and hence we can not use this cri¬ 
terion for the skew-normal distributions. Instead, we use 


where D(lo) is the dispersion function (1M1) and is real at 
u) = rj. See the Appendix [B] for details. The function 
D(rf) can be rewritten as 


D(r]) = 1 + 7T 



Job; A) - job; A) 
(p - r}) 2 


dp, 


(35) 


by performing the integration by parts and remember¬ 
ing job; A) = 0 HU- The Taylor expansion says that 
the numerator of the integrand starts from the quadratic 
term, ( p — rj) 2 , and hence no singularity appears in the 
integrand. A rigorous treatment of the above Penrose 
criterion is found in Ref. H”! • 

The stability criterion (1341) is graphically presented in 
Figl^l The mapped real w axis by D intersects with the 
real D{u> ) axis at u> = rj only, since Im(ZI(aj)) vanishes at 
the unique extreme point. Consequently, we can say that 
the state jo (p; A) is unstable iff the mapped real ui axis 
by D crosses with the negative real axis on the complex 
D(ui) plane. Observing Figj2l the stability threshold of 
the skew-normal distributions, denoted by Ath, must be 
in the interval 1.6 < Ath < 1.7. From symmetry with re¬ 
spect to A, we have another threshold —Ath, and /o(p; A) 
is stable for — A t h < A < A t h- 

The stability threshold can be estimated by precise 
numerical computations. The integral in Eg. (1351) is in an 
infinite interval, and is impossible to perform exactly in 
numerics. To estimate the infinite interval integration, 
we introduce the cut-off P as 


D P (rj) = 1 + 7T 



job; A) - job; A) 


dp, 


(36) 


b - v ) 2 
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FIG. 2. (color online) Nyquist diagrams for the skew-normal 
distributions /o(p;A) with A = 1.5 (green dotted), 1.6 (blue 
dashed) and 1.7 (red solid). Each curve is the mapped real uj 
axis by D, which intersects with the real D(uj) axis at 10 = 17 , 
the unique extreme point. Inside of the curve corresponds to 
the upper half to plane. 



P 

FIG. 3. (color online) Numerical estimation of threshold Ath 
with varying cut-off P (blue circles). The black solid curve is 
the fitting by the least squares method in the interval [10, 100] 
of P, and the red horizontal dashed line is the estimated level 
of Ath = 1.622. 


and observe P-dependence of Ath- Estimated threshold 
with varying P is reported in Figj3l and is fitted by 
1.622 + 1.463/P, where the fitting curve is obtained by 
the least squares method. We hence conclude that the 
threshold is Ath — 1.622 in the limit P —> 00 . 


IV. NUMERICAL TESTS 

We use the semi-Lagrangian code |38] with the time 
slice At = 0.05. The p, space, the ( q,p ) plane, is trun¬ 
cated to (— 7 r, 7r] x [—4,4], and is divided into G x G grid 
points. We call G the grid size. The magnetization is 
zero for the reference homogeneous state /o (p; A), and 
therefore, we simply denote the response magnetization 
as ( M x ,M y ) instead of (Mi jX , 



FIG. 4. (color online) Initial temporal evolutions of M for 
the perturbed initial state f e (q,p; A), (1371) . with e = 10~ 6 and 
A = 1.60,1.61, 1.62, 1.63,1.64 and 1.65 from bottom to top. 
The grid size is G = 512. The vertical axis is in logarithmic 
scale. The stability threshold is in the interval 1.62 < Ath < 
1.63, and is consistent with the estimated value A t h — 1.622. 

It might be worth remarking that the truncation at 
\p\ = 4 does not conflict with the estimation of Ath re¬ 
ported in Fig(3l which requires a larger cut-off. The ref¬ 
erence state fo rapidly decreases as the Gaussian, while 
the integrand in (l36l) slowly decreases as p~ 2 in the large 
\p\ due to existence of the constant fo(rj). 

A. Stability threshold and unstable branch 

The obtained stability threshold is directly examined 
by computing temporal evolution of a perturbed state. 
We prepare the perturbed initial state as 

fe(q,P'A) = /o(p;A)(l + ecosg), (37) 

and use e = ICE 6 . Temporal evolution of M = ( M 2 + 
M 2 ) 1 / 2 is shown in FigjTj and the computed threshold 
Ath is successfully confirmed. 

When the initial state is symmetric with respect to 
p, the nonlinear response theory [lj2j predicts that M 
will be proportional to (A — Ath) 2 hr the unstable branch. 
Numerical simulations captured oscillations of M around 
the predicted levels and the period tends to increase as 
the initial state approaches to the stability threshold [l2| . 
Even the present asymmetric case, the scaling, oscilla¬ 
tions and a similar tendency of periods are observed as 
reported in Figj5] 

B. Linear responses 

We come back to the unperturbed initial distribution 
fo(p; A), and add the external field ([6]). From symmetry 
of the system we set (h x ,h y ) = (h, 0) without loss of 
generality. 
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FIG. 5. (color online) Time averaged M as a function of 
A — Ath for f e (q,p-, A) (13711 with e = 10~ 6 . The time window 
for averages is [1000, 5000]. The grid size is G = 512. The 
green straight line represents M = (A — Ath) 2 /4 for guide 
of eyes. The insets represent temporal evolutions of M for 
the marked points. The horizontal axis represents the scaled 
time t/1000, and the vertical axes 10 5 M and 10 4 M for the 
upper-left and the lower-right insets respectively. 


In order to examine the linear response theory, we set 
h = 10 -5 to be small enough. The normalized responses 
M x /h and M y /h , which are susceptibilities in the limit 
h —y 0, are reported in FigEOfor stable states of A = 1.2 
and 1.6. 

The theoretically predicted levels of responses are in 
good agreements with the numerical experiments in ini¬ 
tial time intervals. The life time of the agreements gets 
longer as the grid size G increases, and is, roughly speak¬ 
ing, proportional to G. We may therefore conclude that 
the theoretically predicted response tensor is valid for a 
long time and that the non-zero off-diagonal response is 
observable if we use a fine grid. 

For the whole stable region of A, the theory is com¬ 
pared with numerical results in Fig(T] We remark that 
the state with A = 0 is the thermal equilibrium state of 
temperature T = 1, and the normalized response M x /h 
coincides with the previously computed Vlasov linear re¬ 
sponse T C /(T — T c ) = 1 plf 0, which is also coincides 
with isothermal linear response (flUl) . We stress that, as 
stated in the end of Sec lII Cl no divergence is observed 
at the stability threshold, which are the left and right 
boundaries of the figure. Another remark is that strength 
of response {M 2 + M 2 ) 1 / 2 /h for A ^ 0 is greater than the 
symmetric case, A = 0. 

One possible explanation for the sign of \ yx is as fol¬ 
lows. We may concentrate for A > 0 without loss of gen¬ 
erality. In this case the negative part of fo(p', A) is larger 
than the positive part around p = 0, and hence the small 
cluster being around p = 0 induced by the external field 
locally has negative total momentum. Consequently, the 
magnetization vector turns to the negative direction of q 
by the external field. 


FIG. 6. (color online) Normalized responses M x /h (left) and 
My/h (right) with h = 10 -5 . A = 1.2 (upper) and 1.6 (lower). 
The grid sizes are G = 128 (green dotted), 256 (blue dashed) 
and 512 (red solid). The black horizontal lines are theoretical 
predictions: M x /h = 2.972 and M y /h = —2.344 for A = 1.2, 
and M x /h = 0.2353 and M y /h = —4.042 for A = 1.6. 



FIG. 7. (color online) Elements of susceptibility tensor as a 
function of the skewness A. Lines are from theory. Points are 
from numerics with the grid size G = 512, and M x and M y 
are computed as averages over the time window [0, 200]. Di¬ 
agonal element Xu (magenta solid/squares), and off-diagonal 
element Xyx (black dashed/circles). The region of A is re¬ 
stricted in the stable interval. 


C. Dependence on external magnetic field 

The present non-diagonalizable susceptibility tensor 
comes from non-zero /q( 0;A), which implies that the 
maximum point p differs from the origin. Thus, we ex¬ 
pect that asymmetric characters of the linear response 
tend to be hidden if the characteristic scale of p-axis, 
width of the separatrix, is larger than the maximum point 
p = p, since the local total momentum in the separatrix 
approaches to zero. 

For the magnetization (M x , M y ) and the external field 

(h, 0), the separatrix reaches to |p| = 2^||M|| + h. The 
magnetization is induced by the external field, and we 
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FIG. 8. (color online) h dependence of susceptibilities for (a) 
A = 1.2 and (b) A = 1.6. Open symbols are for M x /h, and 
filled symbols for M y /h, which are averaged over the time 
windows [0,200] (squares) or [0,100] (circles). The vertical 
black lines represent ht h, and horizontal black lines the linear 
response levels. The horizontal green lines are the zero level. 
The grid size is G = 512. 


have 

11-^11 = h\j (Xxx) 2 + (: Xyx ) 2 - (38) 

Then, we may expect that the asymmetric characters ap¬ 
pear for small h satisfying 

2 

k<hth ’ kth= 4[( Xra ) 2 + (xv*) 2 + 1 ]' (39) 

We report h dependence of susceptibilities in Fig |5] for 
A = 1.2 and 1.6. The normalized responses, M x /h and 
My/h , approaches to the theoretically predicted levels in 
h < /ih, while the off-diagonal response, M y /h , goes to 
zero for larger h. 


V. STATIONARITY AND NONLINEAR 
EFFECTS 


Let us discuss a possible scenario of temporal evolu¬ 
tion with off-diagonal response. First of all, we show the 
fact that the predicted state with non-zero M y is not sta¬ 
tionary by stating that M and h must be parallel in a 
stationary state. 


Jeans theorem 3,0 states that an inhomogeneous 
distribution function is a stationary solution of the 
Vlasov equation if and only if it depends on (q,p) only 
through integrals of the one-particle Hamiltonian system. 
The responded state has non-zero (M x ,M y ) and the in¬ 
tegral is the Hamiltonian 

H = p 2 /2 — M cos(q — a), (40) 

where 

M = {M x + h x ) 2 + (M y + hy ) 2 , tana; = ^ ^ ^ . 

(41) 

Then, for a stationary state f(T-L(q,p)), we have the van¬ 
ishing integral of 


0 = 


sin(q — a) f (Ufa, p))dqdp = M y cosa—M x sina, 


j j [i 

(42) 

since the integrand of the middle term is odd with respect 
to q — a. This equality and the definition of a imply 


My T hy My 

M x + h x M x 


(43) 


and we conclude M and h are parallel. 

As a result, the state predicted by the linear response 
theory is not a stationary state, and hence the system 
does not keep the predicted state as observed in FigEl 
We can point out a similarity of the present phenomenon 
with the nonlinear trapping [s3|. If the Landau damping 
time scale is longer than the so-called trapping time scale, 
then the exponential Landau dam ping stops and a clus¬ 
ter is formed by nonlinear effects [401 ]. In other words, 
the state experiences the linear Landau damping in an 
early time interval, but stops to damp by the nonlinear 
effects. Similarly, the state predicted by the linear re¬ 
sponse theory appears in a short time interval, and then 
disappears. We conjecture that the disappearance conies 
from nonlinearity of the full Vlasov equation. 


VI. SUMMARY AND DISCUSSIONS 

We investigated the response tensor against an asymp¬ 
totically constant external field for spatially homoge¬ 
neous but asymmetric momentum distributions with 
zero-means by using the linear response theory. The the¬ 
ory predicts two interesting characters of the susceptibil¬ 
ity tensor: One is non-diagonalizablility, and the other 
is non-divergence even at the stability threshold. The 
first character implies that the external field added to the 
x-direction induces the magnetization to the y-direction 
even in the simple HMF model. The off-diagonal re¬ 
sponse is not mysterious in our setting, since anisotropy is 
included in asymmetry of momentum distributions. For 
realizing the theoretical setting, we introduced a fam¬ 
ily of skew-normal distributions. After studying stability 






















of the family by the Nyquist method, all the theoreti¬ 
cal consequences are successfully confirmed by direct nu¬ 
merical simulations of the Vlasov equation. We stress 
that the crucial condition for the two characters is non¬ 
zero derivative of the reference state, /q( 0 ) ^ 0 , which 
never happens for symmetric fo(p). One physical exam¬ 
ple of /q(0) ^ 0 can be found in a beam-plasma system, 
whose momentum distribution consists of, for instance, 
a drifting Maxwellian for the beam and a Maxwellian for 
the plasma f2bt . In this example the non-zero derivative 
/o( 0 ) ^ 0 is realized both with and without shifting the 
distribution to set the total momentum zero in general. 
Studying distributions having two or many peaks is a 
future work. 

The state reached by the linear response is neither in 
thermal equilibrium nor in a stationary state, since the 
off-diagonal response is not zero, while the magnetization 
and the external field vectors must be parallel in a sta¬ 
tionary state. The life time of such a state is finite, but 
gets longer as the grid size becomes finer. Thus, we may 
expect that the off-diagonal response can be experimen¬ 
tally observed by using large enough number of particles. 
However, non-stationarity may cause shortness of the life 
time comparing with the symmetric case, and revealing 
the time scale in which the linear response theory is valid 
is remained as another future work. 

Concerning to the above discussion, we remark on va¬ 
lidity of the linear response theory to predict asymptotic 
stationary states. We considered stable reference states, 
and added an external field small enough. Nevertheless, 
the asymptotic stationary states cannot be predicted by 
the linear response theory for asymmetric homogeneous 
initial states. Analogy with the linear Landau damp¬ 
ing might be interesting, which is stopped by nonlinear 
effects. Recently nonlinear equations for magnetization 
moments has been proposed for homogeneous waterbag 
initial distributions in the HMF model under an external 
field [22j] ■ An extension to non-waterbag states possibly 
helps to understand the nonlinear effects and to solve the 
puzzle on the linear response theory. 

In addition to the stable initial states, perturbed un¬ 
stable asymmetric initial states are also studied, and sim¬ 
ilar features are numerically observed with the symmet¬ 
ric case (l 2 | in ordering and oscillations of magnetization 
around the saturated states. Apart from the macroscopic 
variable, looking into difference in distribution functions 
is a remaining work. For instance, the core-halo structure 
f2~il has been observed on the /i space for waterbag initial 
states (U , but it is still unclear if the present asymmetric 
unstable states also yield such structure in the saturated 
states. 

In this article we focused on the asymptotically 
constant external field corresponding to the zero to¬ 
tal momentum, but an oscillating external field of 
cos(wof) (wo £ K.) is also available. Laplace transform 
of the external field provides poles at ui = ±wo, and the 
denominator of susceptibility, |D( 0 )| 2 , is replaced with 
D(±ujo)D(^fuj 0 ) as shown in (IA15I) . where D(oo 0 ) is the 


complex conjugate of D(uiq). As a result, setting wo = V 
where 77 is the maximum point of momentum distribu¬ 
tion, the susceptibility diverges at the stability threshold, 
which satisfies D(rf) = 0. The symmetry is, therefore, not 
essential for the divergence of susceptibility. Even in this 
case, the susceptibility tensor has non-zero off-diagonal 
elements reflecting the asymmetry, see (IA19I) . 
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Appendix A: Derivation of Vlasov linear response 


Let Xq be the Hamiltonian vector field associated with 
the Hamiltonian Hq, 0. which is expressed as 


*0 =Py- 
oq 


(Al) 


Linearizing the Vlasov equation around fo(p), we 
have the formal solution of perturbation fi(q,p,t) as 


h(q,p,t) = - e 


,-(t-a)X 0 


{H^Jojds (A2) 


for the initial condition fi(q,p,t = 0) — 0. The operator 
exp(LAo) acts on a function u(q,p) as 

e tx °u(q,p) = w(^(g,p)), (A3) 


where is the Hamiltonian flow associated with Hq and 
hence Po(q,p) = (q + pt,t) in our setting. We can prove 
the equality 



v(q,p)u{p 0 t (q 1 p))dqdp 



v{ i fl{q^p))u{q 1 p)dqdp 


(A4) 

by changing variables ( q',p') = p 0 t (q,p) and using 
dq'dp' = dqdp from canonical property of p^. Thus, we 
have 


(B) 1 (t) = - If dqdp f B t - s (q,p){H 1 (s),fo}ds, 

J Jli Jo 

(A5) 

where B t (q,p ) = B{p t 0 {q,p)). Performing the Laplace 
transform 03. we obtain 


{B) 1 (uj) = -JJ B u ,(q,p){H 1 (q,ui),fo(p)}dqdp (A 6 ) 

with 


Hi(q, u) = - Mi tX (u) + h x (u) 


COS q 


Mi ,y(u)) + hy[u >) 


(A7) 


Sin q. 












9 


Substituting B = cos q and B = sin q into the linear 
response formula (IA 6 I) , and using the Laplace transforms 
of cos qt = cos (q +pt) and sin q t = sin(g + pt), which are 
respectively 


cos <lu = 7T 


0 -iq 


P l Q 


2 i \p — uj p + ui 


and 


sin q„ 


1 / e~ iq e iq \ 

2 \p — UJ p + U) J 


(AS) 


(A9) 


where uj is the complex conjugate of uj and 

G(u) = [1 - F xx (uj)]F xx (ui) - [F xy (w)] 2 . (A16) 

In particular, the off-diagonal element is written by 


F xy (uj) — 


PV 

7T 2 


f°^-dp — PV 


-oo P- w 


r m dp 

1-00 P + U 




(A17) 


we have the matrix form of 

( = / F xx (u) F xy (uj)\ (/'h x (u))\ 

yMl t y(uj) J \Fy X (uj) Fyy(oj) J \^Ml t y(Uj) J \ky ( W ) / 

(A10) 

The elements of the matrix F are exhibited in (fl9l) . 

To ensure convergence of the Laplace transform urn 
the matrix F(uj) is defined in the upper half uj plane. We 
analytically continue the domain into the whole complex 
uj plane [32|, and the resulting integral is written as 


f /o(p) 

IlPTu 


dp = PV 


/ 


/o(p) 

-oo PF UJ 


dp±S(o;)i7r/o(±w) (All) 


where PV represents the principal value and is the normal 
integral for uj (jL K, and the second term including 

f 0 , Im(w) > 0 

5(w) = \ 1, Im(w) = 0 (A12) 

I 2, Im(w) < 0 


and results to —Im(£)( 0 )) = — tt 2 /o(0) at w = 0 as shown 
in the susceptibility (|25l) . If we consider the oscillating 
external field of cos(wot) (ujq € M), the susceptibility be¬ 
comes 


2X = [1 - F(w 0 )]- 1 F(o;o) + [1 - P 1 (- Wo )]- 1 F(-a; 0 ). 

(A18) 

Thus, for ui o = 77 , where 77 is the unique extreme point 
of f 0 (p), the diagonal elements of susceptibility diverges 
at the stability threshold satisfying D(rj) = 0. Even in 
this case, the oscillating external field gives the non-zero 
off-diagonal element as 


Xxy — 


foi~v )/ 2 


(l + ’"PV / + (n 2 


(A19) 


Appendix B: Nyquist method 


comes from the residues. 

We remark that the linear response (IA 6 I) is rewritten 
as 

CBKH = - ({s u ( g ,p),£i( g ,w)}) o , (A13) 

if we perform integration by parts. The expression (IA13I) 
gives a similar form of the matrix F with the correlation 
matrix C da as 

/ ({cosg w ,cosg }} 0 ({cosg u ,sing }) 0 
^ ^ \({ Sin9 “’ COS<7 })o ({ sinc/ ^- sin(/ }) 

° (A 14 ) 

The matrix F coincides with the correlation matrix C 
as F(ui) = (/3/2)l2 if fo{p) is the Maxwellian with the 
inverse temperature (3. Therefore, the Vlasov linear re¬ 
sponse coincides with the isothermal linear res pon se in 
thermal equilibrium of the homogeneous phase Il3l ll4| . 

In the text we concentrated on response to the external 
field with uj = 0, but a general ui is also available. The 
explicit form of the matrix [1 — F(uj)]~ 1 F(uj) is 


[i-fh]- 1 ^) 


1 ( G(uj) F xy (oj)\ 

D(uj)D(—lj) \~F xy (uj) G(ui) J 

(A15) 


To review the Nyquist method, we restrict ourselves 
in single-peak distributions including the skew-normal 
distributions. Let us define the set R = {D(ui) € 
C | Im(w) > 0 }, where D(ui) is the dispersion function 
(l 23 l) . If this set R includes the origin, then there exists a 
root of the dispersion relation D(u >) on the upper half uj 
plane, and the root corresponds to an exponential grow¬ 
ing mode from the definition of the Laplace transform 

(Ed. 

To study the set R , we investigate the boundary 

dR = {D(ui) S C | Im(w) = 0 }. 

The boundary forms a closed curve, since D{uj) —> 1 
as ui —> ±00. In the limits of uj — > —00 and +00, the 
curve approaches to 1 from the positive and the negative 
imaginary sides respectively, since f$(p) > 0 for p < 77 
and /o(p) < 0 for p > 77, where 77 is the maximum point of 
the single-peak distribution fo(p)- Then, the orientation 
implies that the upper half uj plane is mapped onto the 
inside of the closed curve. The imaginary part of D(u>) is 
proportional to /g(w) for uj real, and vanishes if and only 
if ui coincides with the unique extreme point 77. Thus, 
D(rj) is real and D (77) <0 implies that there is a root of 
D(uj) on the upper half plane (see Figl 3 . 
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